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We use Quantum Monte-Carlo methods to study the ground state phase diagram of a S = 1/2 
honeycomb lattice magnet in which a nearest-neighbor antiferromagnetic exchange J (favoring Neel 
order) competes with two different multi-spin interaction terms: a six-spin interaction Qz that favors 
columnar valence-bond solid (VBS) order, and a four-spin interaction Q 2 that favors staggered VBS 
order. For Q3 ~ Q2 3> J, we establish that the competition between the two different VBS orders 
stabilizes Neel order in a large swathe of the phase diagram even when J is the smallest energy-scale 
in the Hamiltonian. When Qz (Q 2 , J) (Q 2 (Q 3 , J)), this model exhibits at zero temperature 
phase transition from the Neel state to a columnar (staggered) VBS state. We establish that the 
Neel-columnar VBS transition is continuous for all values of Q 2 , and that critical properties along 
the entire phase boundary are well-characterized by critical exponents and amplitudes of the non¬ 
compact CP 1 (NCCP 1 ) theory of deconfined criticality, similar to what is observed on a square 
lattice. However, a surprising three-fold anisotropy of the phase of the VBS order parameter at 
criticality, whose presence was recently noted at the Q 2 = 0 deconfined critical point, is seen to 
persist all along this phase boundary. We use a classical analogy to explore this by studying the 
critical point of a three-dimensional XY model with a four-fold anisotropy field which is known to 
be weakly irrelevant at the three-dimensional XY critical point. In this case, we again find that the 
critical anisotropy appears to saturate to a nonzero value over the range of sizes accessible to our 
simulations. 

PACS numbers: 


I. INTRODUCTION 

Ground states of quantum magnets with S = 1/2 
moments on a two dimensional (2d) bipartite lattice 
(such as square or honeycomb lattices) generally ex¬ 
hibit long-range spin correlations at the Neel wavevec- 
tor qP This T = 0 antiferromagnetic order, en¬ 
coded in a nonzero value of the Neel order parame¬ 
ter vector n, can be destroyed by frustrating further- 
neighborn^ 2 3 or ring-exchange interactions, as well as by 
certain more tractable multi-spin couplings designed^ to 
partially mimic the effect of such frustrating interactions. 
In many examples, the resulting phase has no magnetic 
order and instead exhibits spatial ordering of the bond- 
energy. In such a bond-ordered valence-bond solid (VBS) 
state 1 ^, the singlet projector P(ij) = — Si • Sj + 1/4 of 
two nearest-neighbor spins ( ij) has an expectation value 
that exhibits spatial structure at the VBS ordering wave- 
vector^) K, resulting in a non-zero value for the complex 
VBS order-parameter 'if. 

A standard Landau approach (based on a coarse¬ 
grained free-energy density 11 expressed in terms of pow¬ 
ers of ft and if and their space-time gradients) would pre¬ 
dict that this phase transformation generically proceeds 
either via a direct first-order transition, or via two con¬ 
tinuous transitions separated by an intermediate phase 
which has both orders or no order. Since the latter pos¬ 
sibilities are more exotic, the simplest generic possibility 
within Landau theory is thus a direct first-order tran¬ 
sition. Such first-order behavior is indeed observed in 
square 12 ! an d honeycomb lattice 13 spin models where a 


multi-spin interaction drives the system to a staggered 
VBS state (Fig. 0b)- 

The theory of deconfined quantum critical pointpHUi 
proposed by Senthil et. al. argues that such Landau- 
theory considerations are misleading when the transition 
is towards a state with columnar VBS order (Fig. [iji) 
on the square or honeycomb lattice. Indeed, their 
arguments 14 strongly suggest that such transitions can 
be generically (without fine-tuning any parameter) sec¬ 
ond order in nature. In this alternate approach, one 
writes the partition function as an imaginary-time (r) 
path-integral over space-time configurations n{r,r ), and 
notes that the spatial configuration n(r) on a given time- 
slice admits topological skyrmion textures in spatial di¬ 
mension d = 2. The corresponding total skyrmion num¬ 
ber is conserved during the imaginary-time evolution 
as long as the space-time configuration of n remains 
non-singular. Conversely, when the skyrmion number¬ 
changing operator acts at imaginary time r on pla- 
quette R , it creates a hedgehog defect centered at R, t. 
In this path-integral representation, this hedgehog defect 
carries a Berry-phase 2ttp(R) /q where p(R) =0,1... q— 1 
depends on the sublattice to which R belongs and q = 3 
(q = 4) for the honeycomb (square) lattice case 17 ’®]. 

Remarkably, this phase factor ensures that the trans¬ 
formation properties of T under lattice symmetries are 
identical to those of the complex VBS order parameter 
if for columnar order on both honeycomb and square 
lattice^® 3 ^ The two operators can thus be identi¬ 
fied with each other insofar as their long-distance corre¬ 
lations are concerned (here and henceforth, we refer to 
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ip as the “columnar” order parameter, although ip is also 
non-zero if the system has plaquette VBS order as shown 
in Fig. [T}i for the honeycomb lattice case). The destruc¬ 
tion of Neel order in the ground state can be described 
as a proliferation of such hedgehog defects, providing a 
natural mechanism for a direct transition between Neel 
and columnar VBS orderfSHUD This theoretical descrip¬ 
tion only involves g-fold (g = 3 on the honeycomb lat¬ 
tice and q = 4 on the square lattice) hedgehogs (cor¬ 
responding to T 9 and its Hermitian conjugate), as de¬ 
fects with smaller hedgehog-number carry rapidly oscil¬ 
lating Berry-phases, causing the corresponding terms in 
the action to scale to zero upon coarse-graining. Such 
restrictions on hedgehog charges in space-time configu¬ 
rations of n are best analyzed® in the CP 1 representa¬ 
tion ft = ZaPapzp, where z a is a two-component complex 
field and <x the vector of Pauli matrices. In the CP 1 rep¬ 
resentation, hedgehogs correspond to monopoles in the 
compact U( 1) gauge-field to which the z a are minimally 
coupled 21 23 1 Thus, if the corresponding non-compact 
CP 1 theory (NCCP 1 ) has a second-order transition, and 
if 3-fold (4-fold) monopoles are irrelevant perturbations 
at the corresponding monopole-free fixed point, one ex¬ 
pects that the Neel-columnar VBS transition on the 
honeycomb (square) lattice to be generically continu¬ 
ous, with critical properties in the NCCP 1 universality 
clasJSHUD Conversely, if 3-fold (4-fold) monopoles are 
relevant at the putative NCCP 1 critical point, the sim¬ 
plest scenario is that this leads to runaway flows which 
signal weakly-first order behavior for the Neel-columnar 
VBS transition on the honeycomb (square) lattic^HUD 

To understand the scaling behavior of g-fold monopole 
creation operators in the vicinity of the non-compact CP 1 
critical point, it is instructive to consider a more gen¬ 
eral NCCP^ -1 theory which has TV-component fields z a 
and study the limiting behavior of g-fold monopole per¬ 
turbations in the TV = 1 and TV = oo limits. For in¬ 
stance, four-fold monopoles are known to be irrelevant 
both at TV = il 14 l 19 l 24 tM l and jy = 0 d 14 l 19 l 24 -l, making it 
very likely that they are also irrelevant in the physical 
TV = 2 cas^f^L^. Thus, the Neel-columnar VBS transi¬ 
tion on the square lattice is expected to be generically 
second-order, with critical properties described by the 
NCCP 1 theorjJHHHI. 

The behavior of three-fold monopoles at the noncom¬ 
pact CP 1 critical point is harder to understand from such 
a study of limiting cases. This is because the physical 
TV = 2 case lies between the TV = 1 case where three-fold 
monopoles are rele«ar#- 4 | 19 | 24 | 25 l and lead to a weakly-first 
order transition^, and the TV = o d 14 l 19 l 24 -l limit where 
they are irrelevant. These contrasting behaviors in the 
two limits makes it difficult to argue one way or the other 
concerning the behavior of three-fold monopole perturba¬ 
tions at the TV = 2 critical pointP^HUD A nice summary 
of the expected behavior of the NCCP^ -1 theory with 
g-fold monopoles (including results of numerical simula¬ 
tions) can be found in Ref. [28J 

This theory of deconfined criticality has motivated sev¬ 


eral numerical studies 28 ®^ of model quantum Hamilto¬ 
nians designed 9 Uo host a Neel-VBS columnar transition. 
In parallel work, other studies have tried to access the 
physics of deconfined criticality in three dimensional clas¬ 
sical modelfP® 

On the square lattice (with q = 4), 
QMC simulation j2HHSUM H 3 7l 3 9l 4 ° l 42 l 44 l 45 l find no direct sig¬ 
nature of first-order behavior even at the largest sizes 
studied. This is true both for SU(2) symmetric mod¬ 
els, as well as spin models with enhanced SU(N) sym¬ 
metry, which are expected to exhibit a transition in the 
NCCP^ -1 universality class. Further, critical proper¬ 
ties fit reasonably well to standard scalin g predictions 
for second-order transition j2£ P0 | 33||37 | 39 | 44 | The corre¬ 
sponding values oirjN and rjr), the anomalous exponents 
governing power-law decays of the Neel order parame¬ 
ter n and the VBS order parameter ip, are relatively 
larg d — I— 1^, as expected from the theory of deconfined 
criticality. Additionally, the numerically estimated crit¬ 
ical exponents for large values of TV (using lattice spin 
models with SU(TV) symmetry) approach the limiting val¬ 
ues ob tain ed in a large-TV expansion of the NCCP^' 1 
theor}l® l 39 l 56 l Further, different “designer Hamiltoni¬ 
ans” with different multi-spin coupling d 29 * 34 * 39 ^ yield the 
same estimates for exponents and critical amplitudes. At 
or close to this critical point, histogra ms o f the phase of 
ip exhibit near-perfect U(l) symmetr} l 29 l 33 l 40 ( consistent 
with the idea that the irrelevance of the 4-fold monopole 
insertion operator T 4 implies, via the identification T ~ 
ip, the irrelevance of the 4-fold anisotropy in the phase 
of the VBS order parameter ip. However, in the SU(2) 
case, slow (perhaps logarithmic] drifts with increasing 
linear size L are clearly visible®®' 3 ^® j n certain di¬ 
mensionless quantities which are expected to be scale- 
invariant at a conventional second-order critical point in 
three space-time dimensions — examples include the spin 
stiffness and vacancy-induced spin textures. Since his¬ 
tograms of phase of ip exhibit U(l) symmetry character¬ 
istic of the non-compact theory, it seems plausible that 
these drifts are intrinsic properties of the non-compact 
critical point. This interpretation is supported by the 
fact that Monte-Carlo simulations of a lattice-regularized 
NCCP 1 theorjEHTi] a l so see some drifts that mar other¬ 
wise convincing scaling behavior (it is also possible to 
find different lattice-regularizations that lead to a first- 
order behavioi®l 47 ^. However, at the present juncture, 
there is no detailed understanding of these drifts that 
goes beyond this reasonable guess (see however the recent 
analytical argume nts o f Ref. 1571581) . Finally, we caution 
that some author^®® have also interpreted these drifts 
as either hints of a very weak first order transition, or as 
the signature of a flow towards a new universality class 
different from NCCP 1 . 

What about the honeycomb lattice case (g = 3)? Re¬ 
cent numerical studies of tractable model Hamiltonians 
provide a fairly consistent picture of a direct second- 
order transition between the Neel and the columnar VBS 
statej^- 43 * 44 ! , with numerical estimates of the anomalous 
exponents tjn and pn, correlation length exponent v, and 
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universal scaling function:^! all consistent, within errors, 
with the best estimates for the square-lattice transition. 
Further, slow drifts in spin stiffness analogous to the 
square lattice case, have also been observed at the pu¬ 
tative critical pointful. All this strongly suggests that 
the honeycomb lattice transition is also described by the 
NCCP 1 theory of deconfined criticality. 

However, our recent work has also identified an impor¬ 
tant new feature of the honeycomb lattice transition® 
if the honeycomb lattice transition is indeed described 
by the NCCP 1 theory, 3-fold monopoles must be irrel¬ 
evant at the NCCP 1 critical point. Since ~ ip, this 
would imply that three-fold anisotropy in the phase of 
the VBS order parameter ip is irrelevant at criticality. 
However, it was found® that dimensionless measures of 
this three-fold anisotropy at criticality appear to satu¬ 
rate to a non-zero value as a function of increasing size 
(at least for the sizes at which numerical calculations 
were feasible, which are comparable with those used in 
square lattice studies). The simplest explanation is that 
tripled monopoles are irrelevant with a very small scal¬ 
ing dimension, meaning that the dimensionless critical 
three-fold anisotropy should flow to zero very slowly. If 
one only has access to data over a limited range of sizes, 
it can appear to saturate at a non-zero value. 

The present study aims at clarifying this issue of 
anisotropy, as well as adding some further numerical ev¬ 
idence for the less documented case of deconfined criti¬ 
cality on the honeycomb lattice, relevant for frustrated 
honeycomb lattice spin models 4 '. In this context, we 
note that a recent study® suggests an interesting ex¬ 
perimental realization of deconfined criticality in bilayer 
graphene in magnetic and electric fields, further adding 
to our motivation for studying the Neel-columnar VBS 
transition on the honeycomb lattice. 

We focus here on a numerically tractable model in 
which the nearest-neighbor antiferromagnetic exchange J 
competes with two different multi-spin interaction terms, 
a six-spin interaction Q 3 that favors a columnar VBS 
state (when Q 3 Q 2 ,J), and a four-spin interaction 

Q 2 that favors a staggered VBS (when Q 2 Q 3 , J). 
The deconfined quantum critical point for the model at 
Q 2 = 0 has been studied in our previous worli®, as well 
as in Ref. HU The motivation for perturbing this model 
with the Q 2 term was three-fold: (i) this new energy 
scale (when not too large) will introduce a critical line 
of Q 3 c(Q 2 ^ for the Neel-columnar VBS phase boundary. 
Universality of critical exponents and amplitudes can 
be tested along this critical line® (a) if Q 2 tunes the 
“bare” value of the three-fold anisotropy of the colum¬ 
nar VBS order parameter ip, one could test the behav¬ 
ior of the critical three-fold anisotropy along the phase 
boundary line Q 3 C (Q 2 ); (Hi) the competition between 
the staggered and columnar VBS orders in the regime 
Q 2 ~ Q 3 3> J may reveal exotic physics: the transition 
from staggered VBS order (with maximal winding in the 
valence-bond pattern) to columnar VBS order (with zero¬ 
winding) may proceed through an intervening quantum 



FIG. 1: (Color online) (a) Columnar VBS order on the hon¬ 
eycomb lattice: dark links represent higher values of (Pi) (the 
singlet projection operator on this link) than light links. If 
dark links are instead reinterpreted as representing lower val¬ 
ues of Pi, one obtains a representation of plaquette VBS or¬ 
der at the same wave-vector, (b) Staggered VBS order on the 
honeycomb lattice, where again dark links represent higher 
values of (Pi)- In both figures, we have created different or¬ 
dered domains (represented by different colors and separated 
by dashed lines) by introducing a defect. As already dis- 
cussecP-^, the defect has a spinfull core (a free spin 1/2 sits 
at the domain wall intersection) for columnar/plaquette VBS 
whereas the core is spinless for the staggered VBS. Also shown 
are our conventions for labeling unit cells f, bonds fi belonging 
to unit cells r, and A and B sublattice sites in unit cell r. (c) 
Schematic representations of the 4— and 6 — spins interactions 
terms Q 2 and Q 3 . 


spin-liquid (where no winding sector is favored). 

Before proceeding further, it is useful to summarize the 
key findings of the present work: (i) we establish that the 
transition from Neel to columnar VBS order is continu¬ 
ous for all values of Q 2 , and that critical properties along 
the entire Neel-columnar VBS phase boundary QzciQ?) 
are well-characterized by critical exponents and ampli¬ 
tudes of the NCCP 1 theory of deconfined criticality; (vi) 
the three-fold anisotropy of the phase of the VBS order 
parameter persists all along this phase boundary, with 
slight but perceptible upward drift in its value as Q 2 is 
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increased. To explore the possibility that this may re¬ 
flect the fact that tripled-monopoles are irrelevant with 
a very small scaling dimension, we use a classical anal¬ 
ogy and study the critical point of the 3d XY model 
with a four-fold anisotropy field which is known to be ir¬ 
relevant with a small scaling dimension ^ 5 Sfil. Our results 
for the dimensionless anisotropy on this classical model 
are qualitatively similar to our results for the three-fold 
anisotropy at the Neel- columnar VBS transition: in both 
cases, the anisotropy appears to saturate to a non-zero 
value over the available range of sizes, although, in the 
classical case, one expects it to be irrelevant at the tran¬ 
sition; (in) for Q 3 ~ Q 2 > J, the competition between 
these two different VBS orders does not lead to an inter¬ 
vening spin-liquid phase. Rather, it stabilizes Neel order 
in a large swathe of the phase diagram even when J is 
the smallest energy-scale in the problem. 

The article is organized as follows: in Sec. [TTJ we in¬ 
troduce the J-Q3-Q2 models that we will study and pro¬ 
vide some computational details. In Sec. |III[ we show 
our estimates for the phase boundaries in the (Q 2 ,Q 3 ) 
plane. In Sec. |IV[ we study in greater detail the nature 
of phase-boundary Q 3 c (< 22 ) separating the Neel phase 
and the columnar VBS phase, including the behavior of 
the three-fold anisotropy in the phase of the columnar 
VBS order parameter. In Sec. [Vj we study the classical 
three-dimensional XY model with four-fold anisotropy. 
Finally, we conclude in Sec. VI with a brief discussion 
about possible directions for future work. Some addi¬ 
tional numerical results (on the finite-size scaling anal¬ 
ysis of critical anisotropy, as well as 3d XY model with 
g = 3, 5-fold anisotropic fields) are relegated to Appen¬ 
dices El andlBl 


II. MODEL AND METHODS 

The main focus of our work is the numerical study of 
a model of spin- 1/2 moments on sites of the honeycomb 
lattice, coupled by a nearest neighbor exchange J that 
competes with a four-spin interaction Q 2 and a six-spin 
interaction Q 3 : 


H = 

Hj + Hq 3 4 

-Hq 2 ( 1 ) 

Hj = 

M 



(* 7 > 


Hq 3 = 


P(ij) P(kl) P(mn) “ 1 “ P(jk)P(lm)P(ni) 


(ijklmn) 


Hq 2 = 

-Q 2 ^2 

P(ij)P(lm) P(jk)P(mn) “ 1 “ P(kl)P(ni)i 


( ijklmn ) 


where P(ij) = 1/4 — S,;.S ; is the singlet projector on the 
bond (i,j) and ( ijklmn ) denotes an elementary hexagon 
with vertices labeled cyclically (Fig. 00 - We set J = 1 so 
that all energies are measured in units of J. This model 
is studied using the same techniques as in Ref. 03J for 
both obtaining the ground-state and characterizing its 


physical properties. We summarize them here for com¬ 
pleteness, using the same notations: we use a QMC pro¬ 
jector algorithrrP^ on honeycomb lattices of linear size up 
to L = 60, consisting of L 2 unit cells with two spins cor¬ 
responding to the two-sublattice structure of the honey¬ 
comb lattice. Periodic boundary conditions are imposed. 

Neel order is characterized using the vector order pa¬ 
rameter M = J2r w ith n the local Neel field 
n(r) = Spa — S?b- The unit cell is labeled by f and sub¬ 
scripts A and B refer to the two sites in this unit cell 
located on the different sublattices. The VBS order at 
the columnar wavevector K = (27r/3, — 2n/3) is charac¬ 
terized by the order parameter ip = Tj where V? 

is the local field: 

V r = (P m + e 2ml3 P n + e im / 3 P^)e ,KP , 

with Ppjj (/u = 0 , 1 , 2 ) the singlet projector on one of the 
three bonds g corresponding to the unit cell labeled by 
r (see Fig. [l]). Finally, to quantify the staggered VBS 
order, we follow Ref. □a and use the nematic order pa¬ 
rameter (p = X/r where W? is the local staggered 
VBS order parameter field, written as 

W? = (Pm + e 2vi/3 P ?1 + e 4 ™ /3 P ra ) . 

Note the absence of any r dependent phase factor in this 
definition. This is consistent with the fact that staggered 
VBS order only breaks the symmetry of three-fold rota¬ 
tions, while preserving translational symmetry. 

To detect quantum phase transitions, we consider the 
square of the modulus of the three order parameters of 
interest: ( M 2 ), ( \ip\ 2 ) = (ip^ip), and ( \<p\ 2 ) = For a 

continuous Neel-columnar VBS transition, we expect the 
scaling forms: ( M 2 ) = L~^ 1+r,N ^ f^((Q 3 — Q 3 '^)T 1 / I ' JV ) 
and (\ip\ 2 ) = L-^+^vbs) j^((q 3 _ Q^ c )L 1 / VD ). In writing 
these scaling forms, we assume that the phase bound¬ 
ary is crossed by varying Q 3 at fixed Q 2 and allow for 
two different correlation length exponents vn/d associ¬ 
ated with Neel / columnar VBS correlations at different 
critical values Q ^/ D . We do not quote the scaling form 
for the staggered VBS order as this transition is strongly 
first order. 

We also use the following Binder ratios qm = 
((M 2 ) 2 )/(M 2 ) 2 , p = (|^| 4 )/((|^| 2 )) 2 and g 0 = 

(|£/,| 4 )/ ((|I? 0 | 2 )) to locate the quantum critical points 
where Neel, columnar and staggered VBS orders respec¬ 
tively disappear. The two first Binder ratios are ex¬ 
pected to scale close to a continuous quantum phase 
transitions as qm = 5 m((Q 3 — Q^,)L 1 / ! ' JV ) and g^ = 
9 i>((Qs ~ Q^c)^ 1 ^ 0 ) respectively. Note that both VBS 
Binder ratios are not written in terms of the powers of 
the corresponding VBS order parameter, as this would 
involve computations of 8 —spin correlation functions, for 
which there is no simple expression in the valence-bond 
formalism used in the QMC simulations. In stead , we use 
moments of the Monte-Carlo estimator E ^ 29 l 33 l (respec¬ 
tively E#), whose Monte-Carlo average E (respectively 
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E^) coincides with the quantum-mechanical expectation 
value (ip) ((</>)) of the columnar (resp. staggered) VBS 
order parameter. In all our simulations, we found that 
this correctly reproduces the expected physical behavior 
for moments of ip or <p. 

Close to continuous quantum phase transitions, we 
have fitted our numerical data to the respective scaling 
forms, using polynomial up to second order in most cases 
for the universal functions pM/ip and 

We now introduce the observables related to the phase 
of the columnar VBS order parameter ip. The phase of 
ip distinguishes a fixed columnar (‘Kekule’) pattern of 
bond-energy expectation values from one in which a sub¬ 
lattice of plaquettes hosts a valence-bond resonance (see 
Fig. 0 - Both patterns correspond to a three-fold sym¬ 
metry breaking and lead to order at the same wavevector 
K, but they differ in the phase of the complex VBS order 
parameter ip. In our QMC simulations, we do not have 
access strictly speaking to the phase of ip, but rather 
to the phase 9e^ of the estimator E^ = \E^\ exp (iOe^,)- 
We nevertheless expect that it reflects the behavior of 
the true phase of ip. To address the relevance of 3-fold 
monopole events, we consider the following dimension¬ 
less measure of the anisotropy in the distribution of this 
phase: 


W 3 = J dE^P(E^) cos(36> E J ( 2 ) 

with P(E^) is the normalized probability distribution for 
this quantity as sampled by the Monte-Carlo run. 

It is also p ossible to analyze our data using scaling the- 
orie d 25 ! 26 -^ to capture the finite-size behavior of W 3 near 
criticality. We have used such a scaling analysis to fit our 
numerical data as detailed in Appendix A, but we pre¬ 
fer to display the bare numerical data for the anisotropy 
measure W 3 in Sec. |IVB| in order to avoid any assump¬ 
tion regarding the scaling form obeyed by W3. 


III. PHASE DIAGRAM 


We first present our results on the phase diagram of 
the ground-state of H in the (Q 2 ,Qs) plane. As noted 
earlier, Hj favors Neel ordering, while Hq 2 ( Hq 3 ) favor 
staggered (columnar) VBS order. We can locate two lim¬ 
iting points using results from previous works. The model 
Hj + Hq 3 has been showrPS to host a continuous phase 
transition from the Neel to a columnar VBS state at 
Q 3 c(Q 2 = 0) ~ 1.19. Since Q 2 disfavors columnar VBS 
order, we expect the phase boundary Q 3c (( 32 ) between 
the Neel state and the columnar VBS state to define an 
increasing function of Q 2 , at least for small Q 2 - On the 
other hand, Ref[T51 showed that the model Hj + Hq 2 ex¬ 
hibits a strongly first-order transition from the Neel to 
the staggered VBS state at Q 2 c(Q 3 = 0) ~ 6.4. We ex¬ 
pect the first-order transition to staggered VBS order to 
shift to increasing values of Qo when Q 3 is turned on. 


Neel Binder Ratio 




-2 0 2 4 6 8 10 12 14 16 

Q 2 

Staggered VBS Binder Ratio 



FIG. 2: (Color online) Color maps of Binder cumulants (top 
: Neel Binder cumulant gM, middle: VBS columnar Binder 
cumulant g bottom: staggered Binder cumulant g^,) for dif¬ 
ferent values in the (Q 2 , Q 3 ) parameter space, for a system of 
linear size L = 24. Low values indicate long-range order, while 
high values indicate absence of order. These results allow to 
map the phase diagram where Neel, columnar and staggered 
VBS phases can be identified (green lines are indications of 
approximate phase boundaries). 


A first estimate on the location of these phase bound¬ 
aries is given by the magnitude of the Neel Binder cumu¬ 
lant gM- In our definition of gM, and for a large enough 
system size, a value close to 1 corresponds to a phase 
with antiferromagnetic order, while a value 5/3 corre¬ 
sponds to gaussian fluctuations centered at zero, signal¬ 
ing no magnetic order. At the quantum Neel-columnar 
VBS critical point at Q 2 = 0, the Neel Binder cumulant 
taked^a value ~ 1.42 (which should be universal), lying 
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L = 24 
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FIG. 3: The first-order transition from Neel to staggered VBS 
order is readily identified by the sharps jumps in the corre¬ 
sponding order parameters (up triangles: (M 2 ), down trian¬ 
gles (< i > 2 )), for three values of Q 3 (system size L = 24). 


between these two limiting values. In contrast, close to 
a first-order transition^, this Binder cumulant can take 
values larger than 5/3 on finite-systems. We display the 
magnitude of qm for a system of moderate size L = 24 
in the top panel of Fig. [2j This allows a first estimate 
of the phase boundaries: we clearly observe two transi¬ 
tion lines emerging from the limiting points at Q 3 = 0 
and Q 2 = 0 . The nature of the transitions does not ap¬ 
pear to change, since we observe very high values for qm 
(signaling a first-order transition) for the line emerging 
from Q 2 (Q 3 = 0 ), and intermediate values (between 1 
and 5/3) for the line emerging from Q 3 (Q 2 = 0), sig¬ 
naling a continuous transition. This is confirmed by a 
finite-size scaling analysis in the next section. From this 
study of gM , we also see that antiferromagnetism sur¬ 
vives in the region Q 2 ^Q 3 J. Thus, the competition 
between the two YBS orders does not lead to spin-liquid 
behavior. Rather, it allows antiferromagnetism to set in 
although J is the smallest energy scale in the Hamilto¬ 
nian. The phases where no antiferromagnetism is present 
are naturally expected to host columnar (at low Q 2 ) and 
staggered (low Q 3 ) VBS orders. This is well confirmed 
by the low values (close to 1 ) taken by the columnar and 
staggered VBS Binder cumulants displayed in the middle 
and bottom panels of Fig. [2j 

We now consider more carefully the transition line 
Q 2 c{Qz) between the Neel and staggered VBS order, by 
locating the abrupt first-order jumps in the two order pa¬ 
rameters. An example of these jumps is shown in Fig. [3] 
and the resulting phase boundary is represented as a line 
in Fig. [2j The transition between the Neel and columnar 
VBS transitions deserves a more careful finite-size scal¬ 
ing analysis, which is presented in Sec. m the resulting 
transition line Q 3 c (Q 2 ) is also represented in Fig. [ 2 ] 


IV. NEEL-COLUMNAR VBS TRANSITION 
LINE 

A. Exponents and scaling forms 

We focus here on the nature of the phase-boundary 
between the Neel and the columnar VBS states. Follow¬ 
ing our earlier worlP3 at Q 2 = 0, we locate the point at 
which Neel order is lost using the dimensionless Binder 
ratio 5 m, and the point at which the columnar VBS or¬ 
der turns on using the corresponding Binder ratio g^. 
For four different values of Q 2 , we vary Q 3 to locate the 
quantum phase transition and attempt to collapse the 
Binder ratio data onto the corresponding scaling forms 
(see Sec. 0 - In the analysis, we allow these two scaling 
forms to use different values of Q 3c as well as different 
correlation length exponents vn and vn ■ We also analyze 
the collapse of the modulus squares of order parameters 
(M 2 ) and ( \ip 2 \) according to the forms in Sec. [nj pro¬ 
viding estimates of Q 3c , v N , v D as well as and gp. 



1.46 1.47 1.48 1.49 1.5 1.51 1.52 1.53 1.54 


Qa 



Qa 


FIG. 4: (Color online) Crossing plot of Binder cumulants for 
different system sizes : Neel cumulant <?m (top panel, for 
Q 2 = 0.14) and columnar VBS cumulant g^, (bottom panel, 
for Q 2 = 0.60). Symbols are QMC data, solid lines fits to 
the finite-size scaling form (see text). For the fits, a particu¬ 
lar choice of critical window, minimum system size included, 
and order of universal function has been shown here which 
gave y 2 per degree of freedom equal to 1.53 and 0.97 for the 
plots respectively. For estimates on overall error-bars, refer 
to Tabic Q] 

In Figs. [4] and [5] we provide representative examples 
of the results of such an analysis. Our data all along the 
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(Q 3 -03c)i 1/,/D 


FIG. 5: (Color online) Scaling collapse for different system 
sizes of the Neel order parameter (M 2 ) (top panel, Q 2 = 0.85) 
and columnar VBS order parameter (|V’| 2 ) (bottom panel, 
Q 2 = 20.0) in the critical region. Critical point Q^ c and crit¬ 
ical exponents are obtained by fits to the standard finite-size 
scaling forms (see text). For the fits, again a particular choice 
of critical window, minimum system size included, and order 
of universal function has been shown here which gave y 2 P er 
degree of freedom equal to 1.25 and 1.15 for the plots respec¬ 
tively. For estimates on overall error-bars, refer to Table |T] 


Neel-columnar VBS phase boundary is well-described by 
conventional scaling forms. For ready-reference, we also 
tabulate estimates of the corresponding critical points, 
exponents and amplitudes values obtained using these 
different observables in Table [U 

We find that these estimates of Q?, c at a given value 
of Q 2 agree approximately with each other within sta¬ 
tistical errors. More precisely, the spread in the best-fit 
values of Q^ c obtained from VBS data in two different 
ways (from and (|V ,2 |)) is of the same order as the 
difference in the best-fit Q% c values obtained from scal¬ 
ing collapses of <jy, and gM- The same is true for the 
correlation length exponents and vd at a given value 
of Q 2 . Therefore, we conclude that one can consistently 
account for all the data at a given value of Q 2 in terms of 
a single critical point Q 3 c(Q 2 ) at which Neel order is lost 
and columnar VBS order turns on, with both Neel and 
columnar order parameters controlled by a single correla¬ 
tion length exponent v. Within errors, this estimate of v 
does not exhibit any Q 2 dependence. The anomalous ex¬ 
ponents r \n and i]d are also found to be (^-independent 


within error bars (which are larger for t/d)- Addition¬ 
ally, we note that 77 jv and r\o are close to each other in 
value (although the theory of deconfined criticality does 
not predict that these anomalous dimensions are equal). 
The amplitudes of both VBS and Binder ratios at criti¬ 
cality are also found to be constant within errors along 
the critical line. Finally, we emphasize that all estimates 
of the critical exponents and amplitudes for Q 2 7 ^ 0 agree 
with those found in the case Q 2 = (P. 

Our numerical simulations therefore indicate that the 
entire Neel-columnar VBS transition line belongs to a 
single universality class. Our estimates for the critical 
exponents are very close to the latest estimates for SU(2) 
models on the square latticd^SliaiMl suggesting that both 
honeycomb and square lattice transitions are in the same 
universality class, presumably described by the NCCP 1 
critical theory. This strongly suggests that three-fold 
monopole events are irrelevant at the Neel-columnar VBS 
critical point for a SU(2) model on the honeycomb lat¬ 
tice. 


B. Three-fold anisotropy at criticality 

Given that the entire phase boundary appears to be 
controlled by a single fixed point, it is of interest to in¬ 
vestigate the Q 2 dependence of the three-fold anisotropy 
in the phase of the columnar VBS order parameter if) 
at criticality. To this end, we focus on the histogram 
of measured at and in the close vicinity of our best 
estimate for Q 3 C (Q 2 )- The simplest methodology is one 
that requires the fewest theoretical assumptions about 
the scaling properties of the three-fold anisotropy. In 
this approach, we simply monitor the large-L behavior 
of the dimensionless anisotropy measure W 3 (as defined 
in Sec. 0 for a few values in the vicinity of Q 3 C {Q 2 ) for 
various values of Q 2 . This L dependence is interpreted by 
noting that W 3 tends to zero (respectively to unity) with 
increasing system size deep in the Neel (resp. columnar 
VBS) phase. If three-fold anisotropy is irrelevant at the 
transition, one would expect W 3 to tend to zero for large 
L at the critical point, but increase with increasing L 
when one moves into the VBS phase. 

In Fig. [6j we display the L dependence of this quan¬ 
tity in the vicinity of QzciQz) for three different values 
of Q 2 , two small and one large. From this data, it is 
clear that our earlier finding®] of an apparently non-zero 
large-L limit for this quantity at criticality, remains valid 
all along the Neel-VBS phase boundary, including at the 
largest value of Q 2 studied. This nonzero limiting value 
W 3 C appears to increase slightly with Q 2 , as can already 
be observed in Fig. [ 6 j A critical window around W 3 C can 
be defined by considering the values taken by this di¬ 
mensionless anisotropy in the critical region around Q 3 C 
obtained from the analysis of the previous section. In 
this window, one can attempt a more sophisticated scal¬ 
ing analysis that uses some assumptions about the struc¬ 
ture of the scaling theory for W 3 . This is presented in 







(M 2 ) 

9m 

= (M 4 ) / (M 2 ) 2 

<I0| 2 > 

— 

(\E^)/( 

EA 2 Y 

Q2 

<? 3 c 

Vn 

r\N 

Qzc 

vn 

9m{ 0) 

£? 3 c 

v D 

r/D 

Qic 

VD 

300) 

0.14 

1.496(2) 

0.58(2) 

0.27(3) 

1.496(1) 

0.57(3) 

1.425(2) 

1.483(2) 

0.59(2) 

0.37(3) 

1.491(1) 

0.57(3) 

1.718(5) 

0.60 

2.506(2) 

0.56(2) 

0.31(2) 

2.500(1) 

0.56(2) 

1.427(1) 

2.491(5) 

0.57(3) 

0.23(7) 

2.495(1) 

0.56(2) 

1.721(3) 

0.85 

3.058(2) 

0.55(4) 

0.33(2) 

3.050(2) 

0.56(2) 

1.428(3) 

3.03(1) 

0.60(3) 

0.26(8) 

3.044(2) 

0.56(2) 

1.721(5) 

20.0 

45.3(1) 

0.57(2) 

0.31(3) 

45.27(2) 

0.56(2) 

1.430(2) 

45.0(1) 

0.61(3) 

0.32(6) 

45.18(1) 

0.56(2) 

1.727(1) 


TABLE I: For different values of Q 2 : estimates of critical point, exponent and amplitudes r esultin g from t he finite-size scaling 
analysis of order parameters ( M 2 ), (|'i/’| 2 ) and associated Binder cumulants ( M 4 )/{M 2 ) 2 , (|_E^,| 4 )/(|i?^| 2 ) 2 . Error bars were 
determined from the spread on extracted fit parameters depending on critical window size, minimum system sizes included, 
and degree of polynomial for the universal scaling functions, with y 2 P er degree of freedom always < 1.5. 
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IJL 

FIG. 6 : (Color online) Finite-size dependence of W 3 close to 
the critical point for Qi = 0.14 (top panel), Q 2 = 0.60 (middle 
panel), Q 3 = 20.0 (bottom panel). In each case, we display 
data for one value of Q 3 closer to our estimate of Q 3 C , one 
slightly above and one slightly below. 


Appendix A, and provides independent estimates of W^c 
from fits to a scaling form. These estimates, and the 
resulting conclusions are consistent with those presented 
above from the more direct analysis above. 

We are thus led to two conclusions that appear, at first 
sight, to contradict each other. The first is that critical 
exponents and values of Binder cumulants at criticality 
along the entire phase boundary are compatible with the 
NCCP 1 universality class. The second is that this is ac¬ 
companied by a non-vanishing three-fold anisotropy of 
the phase of -0 at criticality, which furthermore appears 
to vary (albeit slightly) along the critical line. As we 
show in the next section, in the better-understood clas¬ 
sical example of a 3d XY model with weakly-irrelevant 
four-fold anisotropy, the dimensionless anisotropy at crit¬ 
icality again appears to saturate to a non-zero large-L 
limit when studied over a limited range of sizes accessi¬ 
ble to Monte-Carlo simulations. As argued in the next 
section, this suggests a possible rationalization of our 
findings: three-fold anisotropy is indeed irrelevant at the 
Neel-columnar VBS transition, but only very weakly so. 


V. CLASSICAL 3 d - XY MODEL WITH Z 4 
ANISOTROPY ON THE CUBIC LATTICE 

We find it useful to compare this peculiar, apparently 
non-zero large L limit of W3 at criticality to the behav¬ 
ior of an analogous quantity in a much simpler classical 
setting in which one can explicitly tune the bare value 
of the corresponding anisotropy, namely the 3 d — XY 
model with Z 4 anisotropy on the cubic lattice. This 
choice of analogy is dictated by the following consider¬ 
ations: from earlier work, we know that Z 3 anisotropy is 
relevant at the isotropic 3d — XY transition, driving the 
system to a weakly first-order transition, while Z 4 and 
higher anisotropies have all been found to be irrelevant 
at the isotropic XY transition (with Z 4 anisotropy hav¬ 
ing the smallest scaling dimension among the irrelevant 
terms). These conclusions are based on an e-expansion of 
the corresponding field theory 2 ®, Monte Carlo estimates 
of the scaling dimensions of q —fold anisotropy termd^, 
as well as direct numerical simulations of the 3d — XY 
model with Z q > 4 anisotropies (as e.g. in Ref. and of 


Qs = 1.476 
Qs = 1-484 
Qs = 1.492 


Q 2 = 0.14 


Qs = 2.476 
Q 3 = 2.484 
Qs = 2.492 


Qs = 0.6 


Qs = 44.6 
Qs = 45.0 
Qs = 45.4 


Q 2 = 20.0 
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the q = 3 states Potts model 2 - 1 . 

Thus, by adding a Z± anisotropy field /14 to the 
isotropic 3d — XY model and studying the critical point 
as a function of /14, we can study an example of critical 
behavior in the presence of an irrelevant anisotropy which 
scales to zero very slowly (since it has a small scaling di¬ 
mension). This provides us a setting to explore via anal¬ 
ogy the possibility that the nonzero W 3c observed for all 
Q 2 along the Neel-columnar VBS phase boundary could 
reflect the fact that three-fold anisotropy is irrelevant at 
this transition, but has small enough scaling dimension 
that it appears almost marginal (saturating to a non-zero 
value) in the range of sizes accessible to numerics. 

We consider the 3d classical ferromagnetic XY model 
with a anisotropy term, defined by the Hamiltonian 

H = — ^ cos (Off — Opr) — /14 ^ cos( 46 V) ( 3 ) 

(r,r') r 

where (r , r') denotes nearest-neighbor sites on the simple 
cubic lattice and Op are U(l) angular variables £ [0, 2n) 
at site r. This model has a high-temperature paramag¬ 
netic phase where the U( 1) symmetry is unbroken, and a 
low temperature ordered phase where the spins align in 
one of the 4 preferred directions. At h = 0, the model has 
a U{ 1) symmetry which is spontaneously broken in the 
low-temperature phase. To access this physics, we per¬ 
form classical Monte Carlo simulations on simple cubic 
lattice of linear sizes L £ {8,16,24,32,48,64} with pe¬ 
riodic boundary conditions using a combination of local 
Metropolis and Wolff cluster updated. 

We first locate the critical points by a standard scal¬ 
ing analysis for four values of the anisotropy field. To 
this end, we define the vector order parameter m = 
( m x ,m y ) = js J2r( cos ^r),sm(0p)). We measure (|m|) 
(where \m\ = Vm 2 ) and the Binder cumulant B = 
((m 2 ) 2 )/(m 2 ) 2 . We also compute the ratio R of corre¬ 
lation functions at fixed distance R = G p/ 2 /C p/ 4 , where 
Cg. = an< l rt = (£,£,£). The two di¬ 

mensionless observables B and R are expected to satisfy 
the standard scaling forms B = /b((T — T C )L 1 / V ) and 
R = fn((T — T c )L 1 ^) in the vicinity of a second-order 
critical point. Similarly, we also expect the scaling form 

(\m\)=L^f m ((T-T c )L 1 /''). 

We employ this strategy at four values of the 
anisotropy field: /14 = 0,0.5,1.0, 2.0 and present typi¬ 
cal results for these observables in Figs. [7] and [ 8 j Fitting 
to the above forms allows to determine the transition 
temperature T c (h ^) reasonably accurately for each of the 
values of /14 studied. Results of our fits for T c (h±), criti¬ 
cal exponents and amplitudes are given in Table |U| They 
clearly confirm that the universality class of the 3d XY 
model is unchanged by adding a Z\ anisotropic field, i. e. 
it is an irrelevant perturbation at the critical point. Note 
as well how little T c changes as a function of /14. 

Armed with this knowledge, we now study W 4 , a di¬ 
mensionless measure of 4-fold anisotropy in the vicin¬ 
ity of this critical point. We define it analogously to 



T 



T 

FIG. 7: (Color online) 3d XY model with four-fold anisotropic 
field: crossing plot for different system sizes for the Binder 
cumulant B (top panel, for h = 0.05) and correlation ratio R 
(bottom panel crossing plot, for h = 2). 



(t - lyr 1 /" 



FIG. 8: (Color online) 3d XY model with four-fold anisotropic 
field: collapse of the order parameter (|m|) (top panel, for 
h = 0) and the Binder cumulant B (bottom panel, for h = 1), 
according to the scaling forms mentioned in the text. For 
estimates on overall error-bars, refer to Table [II 
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(M> 

Binder ratio B 

Correlation ratio R 

h 

T c 

V 

PIv 

fm( 0) 

Tc 

V 

M0) 

Tc 

V 

M 0 ) 

0.0 

2.201(1) 

0.667(2) 

0.515(1) 

1.106(6) 

2.202(1) 

0.675(10) 

1.2346(5) 

2.202(1) 

0.682(9) 

0.882(2) 

0.5 

2.202(1) 

0.666(3) 

0.51(1) 

1.09(5) 

2.202(1) 

0.676(8) 

1.2365(30) 

2.203(1) 

0.671(1) 

0.882(2) 

1.0 

2.205(1) 

0.665(4) 

0.514(4) 

1.103(10) 

2.204(1) 

0.671(5) 

1.2377(4) 

2.205(1) 

0.67(1) 

0.883(2) 

2.0 

2.212(1) 

0.657(3) 

0.520(5) 

1.13(2) 

2.211(1) 

0.6572(10) 

1.2458(11) 

2.212(1) 

0.665(13) 

0.884(2) 


TABLE II: Estimates of critical temperature, exponent and amplitudes resulting from the finite-size scaling analysis of order 
parameter (|m|), Binder cumulants B = ((m 2 ) 2 )/ (m 2 ) 2 and correlation ratio R = Cl/z/Clu. Error bars were determined 
from the spread on extracted fit parameters depending on critical window size, minimum system sizes included, or degree of 
polynomial for the universal scaling functions, with y 2 P er degree of freedom always < 1.5. 


our definition of W3 for the Neel-VBS transition: W4 = 
f drhP(m)cos(AO m ) with P(fh) the normalized proba¬ 
bility distribution of the order parameter, and 0 rn = 
arctan(m ?/ / m x ) its phase, as measured during the Monte 
Carlo run. In Fig. [ 9 j we show the size dependence of W4 
close to the critical point for two different values of (i 4 
(similar results are obtained for the third non-vanishing 
value of the field studied in our simulations). 

Whereas the anisotropy quantifier W4 increases (to¬ 
wards its limiting value 1) with system size below the crit¬ 
ical temperature, it tends to vanish with system size for 
temperature above T c . At criticality, the anisotropy W4 
appears to be essentially constant (and non-zero), within 



0.02 0.04 0.06 0.08 0.1 0.12 

1/L 

FIG. 9: (Color online) System-size dependence of the 
anisotropy parameter W4 close to criticality for three differ¬ 
ent temperatures: above, below and very close to the critical 
temperature T c . Top panel: h = 0.5, bottom panel: h = 2. 


our range of system sizes for all nonzero /14. We also find 
that this critical value W4 C increases significantly with 
increasing h.4 (see Fig. | 9 ]). A finite-size scaling analysis 
of this behavior, employing some assumptions about the 
finite-size scaling form, is also reported in Appendix [Aj 
and confirms this more elementary analysis. We have 
also studied (see Appendix [B]) the analogous quantities 
for 3 — and 5 —fold anisotropies and find that this unusual 
behavior is specific to the 4 —fold case. 

Our results in the Z4 case for this better-understood 
classical problem are thus entirely analogous to our re¬ 
sults for W3 at the Neel-columnar VBS phase bound¬ 
ary. As in that case, this anisotropy coexists with other 
critical properties being well-fit by standard 3 d — XY 
exponents. Given that Z4 anisotropy is known to be 
weakly irrelevant at the three-dimensional XY transi¬ 
tion, this leads us to suggest that three-fold anisotropy 
is also weakly-irrelevant at the Neel-columnar VBS tran¬ 
sition on the honeycomb lattice. 

VI. OUTLOOK 

We close with a brief discussion of a possible avenue for 
further progress. It would be desirable to have a model 
system where the bare value of the three-fold anisotropy 
in the phase of the VBS order parameter ip could be 
tuned by hand. This would be analogous to tuning /14 in 
the classical three-dimensional XY model. 

To achieve this, we begin with the observation that 
the honeycomb lattice quantum dimer model with ring- 
exchange on hexagonal plaquettes and no inter-dimer in¬ 
teractions is known 65 to order in a plaquette VBS state, 
corresponding to the values ( 2 m+ 1)7 t/3 (m = 0 , 1 , 2 ) for 
the phase of the VBS order parameter ip. The anisotropy 
in the phase of ip in this plaquette-ordered VBS state 
is thus exactly the opposite of the anisotropy in the 
columnar-ordered VBS phase (which corresponds to val¬ 
ues 27 T 77 t /3 for the phase of ip). 

Next, we note that it is possible to write down a six- 
spin interaction term in a SXJ(N) spin model which, for 
large enough N, mimics the ring-exchange term of the 
honeycomb lattice dimer model. This term, given be¬ 
low, is the honeycomb lattice generalization of similar 
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constructions employed recentljBU on the square lattice: 

- R 3 {\(ij)(kl)(mn))((jk)(lm)(ni)\ + h.c.). (4) 

( ijklmn) 

Here, the sum is over all such plaquettes of the hon¬ 
eycomb lattice labelled by ( ijklmn ) with vertices la¬ 
beled cyclically, and \(ij)(kl)(mn)) is the state in which 
(SU (N)) spins i and j form a (SU(1V)) singlet (similarly 
for spins k and l, and m and n). In the large -N limit, 
this reduces to a ring-exchange term on each plaquette. 

With this motivation, we expect that a non-zero R 3 
will counter the columnar phase anisotropy seen at the 
critical point of the SU(2) invariant J — Q 3 model and 
allow us to tune the value of W 3 while leaving other 
critical properties unchanged. Thus, we conjecture that 
the SU(2) invariant J — Q 3 — R 3 model (employing the 
R 3 term defined above) provides a promising setting in 
which one can tune the bare value of the anisotropy in 
the phase of ip, and explicitly check the idea that this 
three-fold anisotropy is a weakly irrelevant variable at 
the Neel-columnar VBS transition. In addition, it may 
even be possible to change the character of the ordered 
state (from columnar to plaquette VBS) if R 3 dominates 
over Q 3 . It should be possible to confirm these ideas us¬ 
ing projector QMC simulations of this J—Q 3 — R 3 model, 
and we hope to return to this in future work. 
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Appendix A: Finite-size scaling analysis of the 
dimensionless anisotropy quantifier 
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FIG. 10: (Color online) Scaling collapse plots according to 
the scaling form W3 = gw 3 ((Q 3 — QscfL 1 ^ 3 ) for the dimen¬ 
sionless anisotropy quantifier W3 in the J — Q2 — Q3 model for 
Q 2 = 0.14 (top panel) and Q 2 = 20 (bottom panel). For the 
fits, similar to Sec. |IV A| a particular choice of critical win¬ 
dow, minimum system size included, and order of universal 
function was taken here which gave y 2 P er degree of freedom 
equal to 1.26 and 1.31 for the plots respectively. 
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Q 2 

Q 3 c 

V3 

9 w 3 (0) 

0.0 

1.183(2) 

0.57(2) 

0.115(6) 

0.14 

1.485(1) 

0.58(1) 

0.120(3) 

0.60 

2.485(1) 

0.56(2) 

0.129(2) 

0.85 

3.027(2) 

0.56(2) 

0.128(3) 

20.0 

45.00(3) 

0.57(1) 

0.134(2) 


TABLE III: Results of finite size scaling analysis for the di¬ 
mensionless anisotropy quantifier W 3 for the J — Q2 — Q 3 
model. Error bars were determined again using the same pro¬ 
tocol as in Sec. IV A of the main text (see Table [I]). 


To supplement the W 3 versus L behavior at fixed Q 2 
that we looked at in the main text, we perform a finite- 
size scaling analysis based on the scaling theory of Lou 
et aP^l. Ref. [2U studied the classical 3d XY model in 
presence of a Z q anisotropy field, which is a danger¬ 
ously irrelevant operator at criticality for q > 4, and 
proposed a scaling form for the dimensionful anisotropy 
order parameter as (m q ) = f mq ((T — T C )L 1 / V ' ! ), 

an extension of the XY order parameter scaling form 
(m) = L~PI V f m {{T — T C )L 1 / 1 '). v q is the exponent asso¬ 
ciated with a length scale below which the order param¬ 
eter distribution appears isotropic, even below T c . We 
have v q > v, as this length scale diverges faster than the 


ferromagnetic correlation length (see the analogy with 
the VBS aniso tropy length scale in the theory of decon- 
fined criticalit y 14 * 1 ^ . Ref. [26] related v q jv to the scaling 
dimension of the anisotropy field, but we note that in a 
recent work this relation was questioned^*. 

J —Q2 — Q3 model — In our case of the dimensionless 
anisotropy order parameter W 3 , we can assume follow¬ 
ing Ref. HD a similar scaling form gw 3 {(Q's — Q 3 c)L lR3 ) 
for fixed Qi, without further assumption on v 3 . Fig. [l0| 
shows examples of this scaling analysis and Tab. Ill sum¬ 
marizes the results of the corresponding fits. 

We see that the critical point Q 3c extracted from the 
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FIG. 11: (Color online) Scaling collapse plots according to the 
scaling form W 4 = < 7 w 4 ((T — for the dimensionless 

anisotropy quantifier W4 in the 3d XY model with 4—fold 
anisotropy field for h = 0.5 (top panel) and h = 2 (bottom 
panel). For estimates on overall error-bars, refer to Table [lV| 


h 

T c 

V4 

fw 4 { 0 ) 

0.5 

2 . 202 ( 2 ) 

0.76(10) 

0.031(4) 

1.0 

2.204(1) 

0.70(2) 

0.062(4) 

2.0 

2 . 211 ( 1 ) 

0.665(20) 

0 . 120 ( 1 ) 


TABLE IV: Results of finite size scaling analysis for 
anisotropy quantifier W4 for the 3d XY model with 4—fold 
anisotropic field. Error bars were determined with the same 
procedure as in the main text (see Table [T|. 

verified (this can be expected as we probably need larger 
systems when anisotropy is stronger). 


Appendix B: 3d XY model with 3— and 5— fold 
anisotropic fields 

Here we show that a nearly-constant critical anisotropy 
is specific to the 3d XY model with 4—fold anisotropic 
field by studying the same model with a 3— and 5—fold 
anisotropy field, replacing the term —I 14 ]TYcos( 46) t -?) by 
— h q cos(qdp) with q = 3, 5 in Eq. [ 3 ] We again com¬ 
pute the Binder cumulant and the anisotropy quantifiers 
W3 and W5 adapting the above definitions. 

q = 3 case — We know that the anisotropy is rel¬ 
evant here, rendering the transition first-order. This is 
clearly seen in the top panels of Fig. [l2| where, for two 


scaling analysis is again in agreement with those gotten 
from other analyses (Sec. IVA). We again find the same 
conclusions as that from visual inspection of W 3 versus 
L behavior: there is a finite value of W^ c = gw 3 ( 0) at 
the critical point for all Q 2 , which furthermore seems to 
slightly increase with Q 2 . Finally, within our precision, 
it is not possible to positively confirm that the extracted 
value of V 3 is larger than v (the two exponents are essen¬ 
tially equal within error bars) : independent of the exact 
relation between the two 2 ^lSil, this indicates that 3—fold 
anisotropy is only very slightly irrelevant, consistent with 
a non-vanishing W?, c within our system size range. 

3 d XY model with f-fold anisotropy field — We per¬ 
form the same analysis for the anisotropy quantifier W4 
of the 3d XY model. In Fig. El we show the scaling 
collapse for W4 with the scaling form IV4 = qwX(T — 
T C )L V" 4 ) as the anisotropy field I 14 is varied. Table 
summarizes the results of the scaling analyses. 


IV 


We find again the critical temperature T c is in agree¬ 
ment with those extracted from other order parameters 
(Sec. 0 and changes very little with / 14 , as already men¬ 
tioned. This analysis confirms that W4 takes a clearly 
non-zero value W 4 C = fw 4 ( 0 ) at the critical point, which 
logically increases with / 14 . In this case, we are able to 
confirm that U 4 > v as found in Ref. [ 2 S 3 except for the 
largest field h = 2 where this relation is only marginally 



2.224 2.226 2.228 2.23 2.232 2.234 2.236 2.238 2.2- 

T 



FIG. 12: (Color online) 3d XY model with 3—fold anisotropic 
field: temperature dependence of the Binder cumulant (top 
panels) and 3—fold anisotropy quantifier W 3 (bottom panels) 
for two different values of /13 = 0.5,1.0. 
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1.9 2 2.1 2.2 2.3 2.4 

T 

FIG. 13: Color online) 3d XY model with 5—fold anisotropic 
field: temperature dependence of the Binder cumulant (top 
panels) and 5—fold anisotropy quantifier W 5 (bottom panels) 
for two different values of h 5 = 1 , 10 . 


different field values, the Binder cumulant show signifi¬ 
cant drifts in the crossing point between two consecutive 
sizes. The bottom panels show the temperature depen¬ 
dence of W3, which also show drifting pseudo-crossing 
points. The clear increase with L of W3 nearest to the 
transition temperature where the pseudo-crossing in the 
Binder cumulant is located indicates that anisotropy is 
relevant at criticality. Note as well how the value of T c 
is substantially modified by /13. 


q = 5 case — Anisotropy is irrelevant here and the 
second order nature of the transition is revealed by the 
nice monotonic crossing behavior of the Binder cumulant 
in the top panels of Fig. |T3] for two different values of h 5 . 
There is no observable drift in T c even when /15 changes 
by a factor of 10 - in fact, one observes that the Binder 
cumulant are essentially the same, indicating the strong 
irrelevancy of 5—fold anisotropy. The bottom panels of 
Fig. [13] show the size and temperature dependence of 
W5, which as expected clearly goes to zero at the critical 
point. We performed a finite-size scaling analys of the 
data (not shown) which yield the expected results, such 
as non-drifting T c , v 5 > v , fw 5 { 0) = 0 and the correct 
3d XY value for v. 
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